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Abstract. We study the connection between the evolutionary replicator dynamics 
, and the number of Nash equilibria in large random bi-matrix games. Using techniques 

p I | of disordered systems theory we compute the statistical properties of both, the fixed 

points of the dynamics and the Nash equilibria. Except for the special case of zero-sum 
games one finds a transition as a function of the so-called co-operation pressure between 
a phase in which there is a unique stable fixed point of the dynamics coinciding with 
a unique Nash equilibrium, and an unstable phase in which there are exponentially 
many Nash equilibria with statistical properties different from the stationary state of 
the replicator equations. Our analytical results are confirmed by numerical simulations 
of the replicator dynamics, and by explicit enumeration of Nash equilibria. 
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1. Introduction 

Replicator equations (RE) describe the evolution of populations of species interacting 
through co-operation and competition. A fitness is assigned to each species, dependent 
on the composition the population and on the availability of resources, and species fitter 
than average increase in concentration while the weight of species less fit than average 
decreases in time. Replicator dynamics (RD) have found widespread applications in 
game theory and economics as well as in population dynamics where an equivalence to 
Lotka-Volterra equations of theoretical biology can be established [1]. 
In the context of evolutionary game theory [2] RE describe games which are played 
repeatedly and where at any time-step the interaction is between individuals randomly 
chosen out of populations of agents. Each player is here taken to follow a pre- 
programmed strategy which they cannot change. Over time a payoff is accrued, and 
agents reproduce as described above, passing on their strategy to their descendents. 
This mechanism can be seen as reinforcement learning from past experience [3]. 
Prior to the launch of this evolutionary approach to game theory the analysis of games 
was mostly concerned with static aspects [4], i.e. with strategically optimal actions by 
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fully rational players and the characterization of potential equilibrium points. It is here 
assumed that the game is played only once. Players choose potentially stochastically 
from a set of available actions, and payoffs are then paid to each participant dependent 
on the decisions of all involved players. The key notion of a Nash equilibrium (NE) 
then refers to a point in strategy space so that no player can increase his payoff by 
unilaterally deviating from this point. 

In this paper we extend existing work on single-population random replicator dynamics 
[5] in which a transition between stable and unstable regimes has been found as a 
function of a so-called 'co-operation pressure', and study the case of multi-population 
models, i.e. interaction between distinct populations of species. These correspond to 
so-called 'asymmetric' games [1] in which there is more than one type of player, such as 
in the game known as 'the battle of the sexes' in which male and female types of agents 
have different strategy sets at their disposal. In particular we study the stability of two- 
population replicator systems, and the relation to the number of NE in the corresponding 
bi-matrix games. We here extend the work of Berg et al [6], who computed the number 
and statistics of NE of matrix games with random payoff matrices in the absence of co- 
operation pressure and who identified an exponentially large number of NE for all types 
of payoff correlations considered, except for the case of zero-sum games in which there 
is a unique NE. Our work thus extends the statistical mechanics approach to replicator 
systems to the case of multiple populations, and establishes a connection between the 
earlier works of [5] and [6]. 

2. Model definitions 

We will consider bi-matrix games between two players, one of type X and one of type 
Y, with strategy sets T, x and Ti y , and payoff functions fj, x ,fi y : Sa; x £ y — > KL We will 
here restrict to the case in which |E X | = \H y \ = N. Strategies available to players of 
type X are labelled by % G {1, . . . , N}, the ones available to F by j G {1, . . . , N}. 
The extension to more general cases with different numbers of available strategies 
is straightforward. Mixed strategies are then probability distributions over H x and 
T, y respectively, described by vectors x = (xi, . . . ,xn) and y = (y 1 , . . . , y N ), with 
< Xi < 1 and X^ 2 -* = 1> an d similarly for the {yj}- Pure strategies are recovered as 
unit vectors. If player Y plays mixed strategy y then the expected payoff for X's pure 
strategy % reads vf(y) = J2j t i (hj)llj- If both players play mixed strategies the expected 
payoff for X is denoted by i/ x (x, y) = J2i x i u i(y) an d similarly for v y (x., y). A NE is 
then a point (x*, y*) so that v x (x*, y*) = max x z/ :r (x, y*) and ^(x*, y*) = max y i/ 3 .(x*, y), 
and may be characterised by the following conditions 



The inequalities here ensure that no pure strategy delivers a higher payoff than 
obtainable at the Nash point, and the equalities enforce that all pure strategies which 



<[^(y*)-^(x*,y*)]=0, ^(y*)-^(x*,y*)<0, V* 
y) [if (x*) - i/„(x*, y*)] = 0, z/jV) - i/„(x*, y*) < 0, Vj. 



(1) 
(2) 



Two-population replicator dynamics and number of Nash equilibria in matrix games 3 



are actually played with non-zero probability yield the same payoff, and that all other 
pure strategies are not employed. 

RE of evolutionary game theory assume large populations of X and F-type players, 
respectively, with each individual playing a pre-specified pure strategy. The replicators 
are hence the pure strategies of the game under consideration and are copied without 
error from parent to child. Xi(t) denotes the relative concentration of pure strategy i 
in the X-population, and similarly for Hj(t). The replicator equations describing the 
evolution of the X and Y populations are then as follows 

±i(t) = Xi (t) [*f (y(t)) - K x (t)] , in(t) = yi (t) [**(x(t)) - K y {t)] (3) 

where the Xi(t) and yj(t) are now time-dependent variables, and where K x (t) = 
Y^i x i('t) l 'i(y('t))j an d similarly for K y (t) are the average payoffs of the X and Y 
type players respectively. The replicator equations preserve the overall normalisation 
Yli x i(t) = YljVj(t) = 1 in time. In terms of population dynamics the weights {xi(t)} 
and {yj(t)} correspond to relative concentrations of species, and the payoff functions vf 
an v v - may be seen as their respective fitnesses. In a biological setting RD thus describe 
the temporal evolution of populations of species, where the concentrations of species 
fitter than average grow in and where all other species decrease in relative numbers. In 
the remainder of the paper we will use both the game theoretical and the population 
dynamical language synonymously. Note that in the two-population system species of 
type X interact only with species of type Y and vice versa. 

We will in the following be concerned with replicator systems where all entries of the 
payoff matrices are Gaussian random variables. Following Peschel and Mende [1] as well 
as [5] we will also introduce so-called 'co-operation pressures' u x > and u y > 0, and 
will consider payoff functions of the form 

( x > y) = -2u x Xi + a iiVii Vjfc y) = - 2u yVj + Yl b ^ Xi ' ^ 

3 i 

The role of u x , u y will be clarified below. For u x = u y = one recovers the cases studied 
in [6] (with payoff matrices = (J, x (i,j), b^ = fj, y (i,j)). Both the static and dynamical 
properties of bi-matrix games are invariant under global rescaling and shifts of all payoff 
matrix elements. Without loss of generality we may therefore assume that the Gaussian 
variables and bjj are drawn from the following statistics 
— — 1 r 

a % = b % = ]y> a ^ b kl = dirfjkjj (5) 

(with 7TT an average over the distribution of payoffs). The scaling with iV is here 
introduced to guarantee a well defined thermodynamic limit. The parameter T 
characterises the correlations between the payoff matrices for the two different types 
of players. If F = —1 then = —bji so that the resulting game is a zero-sum game 
(at vanishing co-operation pressures), corresponding to a prey-predator relation in the 
population dynamical setting. For T = one has uncorrelated payoff matrices, and in 
the fully symmetric case T = 1 the two interacting players always receive equal payoff. 
Note also that for convenience we will re-scale the Xi and y 3 - such that they obey the 
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normalisations ^2 i=1 Xi(t) = ^2j = iUj(t) = N at all times t. The co-operation pressures 
u x and u y finally control the growth of individual species, large values of u x drive the 
configuration x into the the interior of the simplex defined by Yli %i = N,0 < Xi < N, 
and similarly for u y , see the book by Peschel and Mende for further details [1]. Indeed, 
inspection of Eqs. (3) for u x ,u y — > oo shows that the RD leads to the fixed point 
Xi = yj = 1 for all species when started from non-zero concentrations, corresponding 
to full co-operation and maximal diversity. In a game theoretic setting large values of 
u x ,u y favour the use of mixed strategies as opposed to pure strategies located at the 
corners of the above simplices. 



3. Dynamics: generating functional analysis 

On here directly addresses the dynamics described by the replicator equations, and 
formulates an effective theory for macroscopic dynamical order parameters. The starting 
point of the analysis is the moment generating functional, defined as 

Z[^»] = ^exp^ J *{f]^(^i(0 + E^(*)l/i(*)}^, (6) 

where ((. . .)) denotes an average over trajectories of the system, i.e. over solutions of 
the replicator equations (3). As usual for disordered systems a closed, but implicit 
set of equations describing the temporal evolution of a small number of disorder- 
averaged macroscopic order parameters can be derived. In the thermodynamic limit 
these observables turn out to be given by the correlation functions C x (t,t'),C y (t,t'), 
the response functions G x (t,t'),G y (t,t') and the Lagrange parameters K x (t),K y (t). The 
correlation and response matrices for the X-population are given by 

c x (t,f) = jim N-^T&M^m, G x (t 1 if) = ^N- 1 ^(^^yr) 

and analogous definitions for C y , G y apply. These order parameters are to be determined 
self-consistently as averages C x (t,t') = (x(t)x(t'))^ G x (t,t') = (5x(t)/5K x (t))+ (and 
analogously for C y , G y ) over realisations of the following pair of coupled stochastic 
effective processes 

x(t) = - x(t) 2u x x(t) - T f dt' G y {t,t')x(t') - K x (t) + r]x{t) (8) 
y(t) = - y(t) hu y y{t) -Y [ dt' G x (t, t')y(t') - n y {t) + Vy (t)] (9) 

L Jt Q J 

with t the starting point of the dynamics. The notation (. . )^ refers to an average over 
the effective process, i.e. over realisations of the noise variables {f] x {t)} and {%{t)}. The 
covariances of these noise variables are given by (r) x (t)r) x (t') ^ = C y (t,t'), {Vy(t)r}y(t') ^ = 
C x (t,t') with no correlations between r] x and r] y . Finally, the Lagrange multipliers 
{ K x(t)} and {tt y (t)} have to be chosen such that the constraints (x(t))^ = 1 and 
= 1 are fulfilled at any time t >t . 
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Further progress can be made by assuming a fixed point of the replicator equations, i.e. 
by inspecting time-independent solutions x(t) = x, y{t) = y of the effective processes 
(leading to constant correlation functions C x (t,t') = q x ,C y (t,t') = q y and to stationary 
response functions G x (t), G y {r)). Similarly to [5] one derives the following equations 
characterising such fixed point states 

Xx(2u x - Yx v ) = go(Ax) Xy(2u y - Yxx) = go(\) 

q y 1/2 (2u x - T Xy ) = gi(A x ) q x l/2 (2u y - V Xx ) = gi(\) 

(q x /q y )(2u x - Txyf = ^ 2 (A X ) (q y /q x )(2u y - Txxf = #2(A y ) (10) 

with A x = k x / ' jq- y , A y = Ky/y/ch and g n (A) = dz^£ L (A - z) n for n E {0, 1, 2}. 
Xx is here the integrated response Xx = f dtG x (r), and similarly for Xy A linear 
stability analysis shows that such fixed points become unstable at the point at which 
(XxXy) 2 — leading to an unstable and non-ergodic phase at low co-operation 

pressures. For u x = u y , which we will mostly consider in the following, one finds 
Xx = Xy an d q x = Qy, and Eqs. (10) as well as the stability condition reduce to those of 
a single-population replicator system studied in [5]. In the unstable phase solutions of 
(10) can no longer be expected to describe the stationary states of the RE accurately. 



4. Statics: replica analysis 



a ij = bji (r = 1) the replicator equations (3) 
Xi (d Xi 7i(x, y) — n x (t)) and similarly for ijj, with 



The starting point of the replica analysis of the statics of the model is the observation 
that for symmetric couplings, i.e. 
can be written in the form Xi = 

W(x,y) = 1X^=1 [xiVjioij + bp)] - u x J2" =1 xj - u y Y^ =1 y 2 j, so that the stationary 
states of the RD at Y — 1 correspond to extrema of H. The computation of these is 
straightforward and based on the evaluation of —f3f = hin^^oo iV -1 ln Z([3). f is here 
the disorder- averaged free energy density at temperature T = Z(f3) stands for the 
partition function corresponding to 7i: 



Z(P) 



N 



N 



H I dx U 



J'=l 



dy 3 



N 



N 



N 



N 



-/9W(x,y) 



(11) 



We will not report the details of the straightforward replica calculation, but will only 
note that a replica-symmetric ansatz leads to a set of equations identical to (10) at zero 
temperature and in the thermodynamic limit (with T = 1). For T < 1 no Lyapunov 
function Ti can be found and the replica approach is inapplicable. 



5. Annealed calculation of the number of Nash equilibria 

Finally, the number and statistics of the NE of the corresponding bi-matrix game may 
be computed by direct integration over phase space enforcing conditions (1,2) through 
suitable delta- and step-functions. It is here convenient to set = Xi if Xi > and 
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Xi = —2u x Xi + dijUj — k x if Xi = 0, and similarly for yj [6]. The above conditions (1) 
may then be written as 



Jf(x,y) = XiB(-Xi) - ^-2uxiQ(xi) + ^ a ij y j e(y j ) - K^j =0 



(12) 



and (2) translates into Jj(x, y) = with an analogous expression ij(x, y). ©(•) is the 
step-function. The number of NE at payoffs k x and n y is then given by 

A/\k x , K y ) = j DSiDyS ^x&ixi) - N^j 5 ^yMVi) - ivj 

*I[6(I?(%y))I[6(I?(3c,y)) \detD\. (13) 

i j 

det-D is a normalising determinant. Performing the disorder- average in an annealed 
approximation, one converts the computation into a saddle-point problem in the 
thermodynamic limit. We set for simplicity and find exponential 

domination at equal payoff and limjv^oo ^\o.H{k) = S(k), where 

S(k) = extr {£ R q S^E - TR 2 + 2qq + 200 + 2 In H (--^ 



V^qJ y/(2u - TR) 2 + 2qq 



, (2u-YR)K/q + E \ ( k 2 (E+(2u-TR)K/q)< 

xH = x exp -— + 



^/q-\2u-TR) 2 + 2q) * \ 2q 4q + 2(2u - TR) 2 /q 

We have here set abbreviated H(x) — \ (l — erf (x/ V2)) , and g denotes the contribution 
from the normalising determinant, which we do not report explicitly. For u = we 
recover the result of [6]. Extremisation of (14) leads to an annealed upper bound of the 
logarithmic number of NE at payoff k. 



6. Results 



Our results are summarised by figures 1 and 2. Fig. la shows the phase diagram of 
the bi-population replicator model at different symmetry parameters T. For T > — 1 
one finds a transition line, with unique fixed points of the RE at large co-operation 
pressures to the top-right, and an unstable regime below. For T = — 1 no transition is 
observed, in this case there is one unique fixed-point of the dynamics for any (u x ,u y ). 
In the remaining figures we restrict to the case u x = u y = u. The transition between the 
stable and unstable phases then occurs at co-operation pressure u c (T) = (1 + r)/(2y / 2). 
In Fig. lb we depict typical curves S(k) of the entropy of NE at payoff k as obtained 
from Eq. (14). These curves typically show a maximum S m (u, T) at intermediate values 
of the payoff k, indicating the number of dominating NE. If S m > a large number 
of Nash equilibria is present in the thermodynamic limit, while for S m < NE are 
exponentially suppressed. If S m — a single NE prevails [6]. The inset of Fig. lb shows 
a comparison between the analytical results for S(k) at a fixed combination of u and V 
with numerical data from an explicit enumeration of NE. We here use a direct support 
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Figure 1. 

Left: (Colour on-line) Phase diagram of the two-population replicator system. Lines 
separate the stable phase (up-right) from the instable one, and are plotted for 
r = 1, 0.5, 0, —0.5 from top to bottom. 

Right: (Colour on-line) Main panel: entropy of number of Nash equilibria for the 
uncorrelated game r = at u = 0, 0.1, 0.2, 0.3 (from top to bottom at the maximum). 
Inset: comparison with numerical data from support enumeration (it = 0.1). Upper 
curve is the theoretical line, others are from support enumeration of single-population 
models at N = 14, 16, 18, 20 (bottom to top). A running average has been performed 
to smoothen numerical data. 

enumeration scheme, which we apply to an equivalent single-population game in order 
to reduce the required computing time. Accessible system sizes are limited to iV rs 20, 
and while results are consistent with the theory, finite size effects can be significant. 
Fig. 2a shows the entropy S m of the dominating NE as a function of u at fixed values 
of the symmetry parameter T. One finds a large number of NE (S m > 0) for u < u c (r), 
and a unique NE above u c (S m = 0). At u c the spectra S(k) take their maxima at zero 
payoff k = 0. The transition point u c (T) coincides with the one separating the stable and 
unstable phases of the RD. Results from an explicit enumeration of NE are consistent 
with this transition although direct quantitative comparison between simulations and 
theory appears difficult due to finite size effects. An example of raw data obtained from 
enumeration of NE is shown in the inset of Fig. 2a. 

The saddle point extremisation of (14) allows one also to compute the statistics of the 
dominating NE (in the annealed approximation) and to compare with the stationary 
states of the RE. We here focus on the order parameter g -1 which serves as a measure 
of the diversity of the eco-system [7], with large values of q~ l corresponding to many 
surviving species (equivalently many pure strategies played with non-zero probability). 
Measurements of g -1 in the stationary states of the RE and numerical results for the 
NE are shown in Fig. 2b. The latter are here obtained using a repeated Lemke- 
Howson algorithm [8]. Note that the fixed-point theory of the RE applies only above 
u c but has been continued below as dashed lines. Quantitative deviations between 
numerical and theoretical results for NE are due to finite-size or sampling effects and 
the annealed approximation of the theory. The analytical and numerical results show 
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u u 
Figure 2. 

Left: (Colour on-line) Entropy S m of Nash equilibria as a function of u, for T = 
(left curve) and T = 1 (right curve). Solid lines from theory, symbols from explicit 
enumeration of single-population systems with N = 20, averaged over 1000 samples. 
Inset: Raw data for S(k) at T = 1 for u = 0.9, 0.8, . . . , 0.0 from left to right. 
Right: (Colour on-line) Diversity parameter q" 1 at T = —1,0, 1 (from left to right). 
Solutions of (10) depicted as solid lines where replicator equations have stable fixed 
points, continued into the instable phase as dashed lines. Open symbols are results 
from simulations of bi-replicator systems with 2N = 1000 species. Dashed-dotted 
curves show theory for Nash equilibria for u < u c (T), open symbols are numerical 
results for typical NE as obtained from iterated Lemke-Howson algorithm applied to 
bi-matrix systems of size 2N = 100, averaged over 200 samples. 

that the diversity of species in the NE coincides with that at the fixed points of the 
RD above u c (T), but that they differ from each other below the transition. Very similar 
results are found for the number of surviving species. We note that conditions (1,2) 
are necessary but not sufficient for the stability of fixed points of the RE so that stable 
fixed points are always NE but not necessarily vice versa. Results thus indicate that NE 
are statistically distinct from potential stable attractors of the dynamics below u c (and 
more numerous). Note also that below u c marginally stable fixed points are suppressed 
for T < 1 [5] and hence volatile, possibly chaotic behaviour is observed in simulations. 
For T = 1 the dynamics converges to a fixed point also below u c , but is sensitive to 
initial conditions (see also [5, 9]). 

7. Summary and concluding remarks 

Our results may be summarised as follows: for games different from the zero-sum type 
(i.e. for any T > —1) one finds a dynamical instability of the fixed points of the two- 
population replicator equations at u c (T) = (1 + r)/(2y / 2). For T = 1 this instability 
coincides with an AT-instability of the replica symmetric solution of the statics [10]. 
Above u c all three approaches (dynamics, NE and replica theory where applicable) 
lead to the same order parameters, describing an ergodic state, in which there is one 
unique NE which coincides with the unique fixed point of the RD (and with the unique 
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extremum of 7i for T = 1). Below u c ergodicity is broken, the stationary state of the RD 
is no longer unique, and there is an exponential number of NE. The statistics of the NE 
differs from those of the stationary states of the dynamics. The observations generalise 
the results of [6] , which were concerned with the case of vanishing co-operation pressure 
u. At u = the RD is instable for all T > — 1, and marginally stable for zero-sum 
games, T = — 1. Except for the latter case the system is therefore in the non-ergodic 
phase of the RE, in line with the reported exponential number of NE in [6]. At T = — 1 
one is precisely at the phase transition at u c (T = —1) = and finds a zero-entropy of 
dominating NE. 

In conclusion we have investigated bi-matrix games and two-population replicator 
systems with tools of disordered systems theory, and have have computed the number 
and typical properties of NE and of the fixed points of the corresponding replicator 
equations. Several extensions of the present work can be considered, including cases of 
correlations between different rows and/or columns of the payoff matrices, populations 
with different numbers of strategies at hand and systems with more than two types of 
players, as well as the study of other stability concepts of evolutionary game theory in 
the context of random replicator systems. 
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